In [1]:

    

import sympy
import sympy.physics.mechanics as mech
import scipy.optimize

sympy.init_printing()
sympy.physics.vector.init_vprinting()
%matplotlib inline


    

In [2]:

    

alpha, beta, eta, P, Q, R, r, alpha_l, beta_l, mu = \
    sympy.symbols('alpha, beta, eta, P, Q, R, r, alpha_l, beta_l, mu', real=True, imaginary=False)
V_T, s, rho, c, c_0 = sympy.symbols('V_T, s, rho, c, c_0', postive=True, real=True, imaginary=False)
C_L_bar = sympy.Symbol(r'\bar{C_L}', real=True, imaginary=False)
C_D_bar = sympy.Symbol(r'\bar{C_D}', real=True, imaginary=False)
t = sympy.Symbol('t')

sub_adv = {V_T(t) : mu*s*eta(t).diff(t)}
mu_sub = {mu:V_T(t)/(s*eta(t).diff(t))}

#inertial frame
frame_i = mech.ReferenceFrame('i')

# body frame
frame_b = mech.ReferenceFrame('b')
frame_b.set_ang_vel(frame_i, 
    P(t)*frame_b.x + Q(t)*frame_b.y + R(t)*frame_b.z)

# wind frame
frame_w = frame_b.orientnew(
    'w', 'Body', (-alpha(t), beta(t), 0), '231')

# wing frame
frame_l = frame_b.orientnew(
    'l', 'Axis', (eta(t), frame_b.x))

# local wind frame
frame_m = frame_l.orientnew(
    'm', 'Body', (-alpha_l(t), beta_l(t), 0), '231')

# center of mass
point_cm = mech.Point('cm')
point_cm.set_vel(frame_b, 0)
point_cm.set_vel(frame_i, V_T(t)*frame_w.x)

# point on wing
point_p = point_cm.locatenew('p', r*s*frame_l.y)
point_p.set_vel(frame_l, 0)
point_p.v2pt_theory(point_cm, frame_i, frame_l)


    

    
Out[2]:





$$V_{T}\mathbf{\hat{w}_x} -  r s \left(- Q \operatorname{sin}\left(\eta\right) + R \operatorname{cos}\left(\eta\right)\right)\mathbf{\hat{l}_x} + r s \left(P + \dot{\eta}\right)\mathbf{\hat{l}_z}$$

In [3]:

    

alpha_l_expr = sympy.atan(point_p.vel(frame_i).dot(frame_l.z)/
    point_p.vel(frame_i).dot(frame_l.x)).simplify()
alpha_l_expr


    

    
Out[3]:





$$\operatorname{atan}\left(\frac{r s \left(P + \dot{\eta}\right) + \left(\operatorname{sin}\left(\alpha\right) \operatorname{cos}\left(\beta\right) \operatorname{cos}\left(\eta\right) - \operatorname{sin}\left(\beta\right) \operatorname{sin}\left(\eta\right)\right) V_{T}}{r s \left(Q \operatorname{sin}\left(\eta\right) - R \operatorname{cos}\left(\eta\right)\right) + V_{T} \operatorname{cos}\left(\alpha\right) \operatorname{cos}\left(\beta\right)}\right)$$

In [4]:

    

beta_l_expr = sympy.asin(point_p.vel(frame_i).dot(frame_l.y)/
    point_p.vel(frame_i).magnitude())


    

In [5]:

    

#beta_l_expr


    

In [6]:

    

q_expr = (rho*point_p.vel(frame_i).magnitude()**2/2).collect(r).simplify()
q_expr


    

    
Out[6]:





$$\frac{\rho}{2} \left(r^{2} s^{2} \left(\left(Q \operatorname{sin}\left(\eta\right) - R \operatorname{cos}\left(\eta\right)\right)^{2} + \left(P + \dot{\eta}\right)^{2}\right) + 2 r s \left(\left(Q \operatorname{sin}\left(\eta\right) - R \operatorname{cos}\left(\eta\right)\right) \operatorname{cos}\left(\alpha\right) \operatorname{cos}\left(\beta\right) + \left(\operatorname{sin}\left(\alpha\right) \operatorname{cos}\left(\beta\right) \operatorname{cos}\left(\eta\right) - \operatorname{sin}\left(\beta\right) \operatorname{sin}\left(\eta\right)\right) \left(P + \dot{\eta}\right)\right) V_{T} + V^{2}_{T}\right)$$

In [7]:

    

Lp_expr = -C_L_bar*rho*q_expr*c(r)/2*frame_m.z
Lp_expr


    

    
Out[7]:





$$-  \frac{\bar{C_{L}} \rho^{2}}{4} \left(r^{2} s^{2} \left(\left(Q \operatorname{sin}\left(\eta\right) - R \operatorname{cos}\left(\eta\right)\right)^{2} + \left(P + \dot{\eta}\right)^{2}\right) + 2 r s \left(\left(Q \operatorname{sin}\left(\eta\right) - R \operatorname{cos}\left(\eta\right)\right) \operatorname{cos}\left(\alpha\right) \operatorname{cos}\left(\beta\right) + \left(\operatorname{sin}\left(\alpha\right) \operatorname{cos}\left(\beta\right) \operatorname{cos}\left(\eta\right) - \operatorname{sin}\left(\beta\right) \operatorname{sin}\left(\eta\right)\right) \left(P + \dot{\eta}\right)\right) V_{T} + V^{2}_{T}\right) \operatorname{c}\left(r\right)\mathbf{\hat{m}_z}$$

In [8]:

    

Dp_expr = -C_D_bar*rho*q_expr*c(r)/2*frame_m.x
Dp_expr


    

    
Out[8]:





$$-  \frac{\bar{C_{D}} \rho^{2}}{4} \left(r^{2} s^{2} \left(\left(Q \operatorname{sin}\left(\eta\right) - R \operatorname{cos}\left(\eta\right)\right)^{2} + \left(P + \dot{\eta}\right)^{2}\right) + 2 r s \left(\left(Q \operatorname{sin}\left(\eta\right) - R \operatorname{cos}\left(\eta\right)\right) \operatorname{cos}\left(\alpha\right) \operatorname{cos}\left(\beta\right) + \left(\operatorname{sin}\left(\alpha\right) \operatorname{cos}\left(\beta\right) \operatorname{cos}\left(\eta\right) - \operatorname{sin}\left(\beta\right) \operatorname{sin}\left(\eta\right)\right) \left(P + \dot{\eta}\right)\right) V_{T} + V^{2}_{T}\right) \operatorname{c}\left(r\right)\mathbf{\hat{m}_x}$$

In [9]:

    

assump_subs = {P(t):0, Q(t):0, R(t):0, alpha(t):0, beta(t):0, c(r):c_0}
assump_subs.update({
    beta_l(t): beta_l_expr.subs(assump_subs),
    alpha_l(t): alpha_l_expr.subs(assump_subs)})


    

In [10]:

    

assump_subs


    

    
Out[10]:





$$\left \{ P : 0, \quad Q : 0, \quad R : 0, \quad \alpha : 0, \quad \alpha_{l} : \operatorname{atan}\left(\frac{r s}{V_{T}} \dot{\eta}\right), \quad \beta : 0, \quad \beta_{l} : 0, \quad \operatorname{c}\left(r\right) : c_{0}\right \}$$

In [11]:

    

Lp_assump = Lp_expr.subs(assump_subs).to_matrix(frame_l).subs(assump_subs).subs(sub_adv).simplify()
Lp_assump


    

    
Out[11]:





$$\left[\begin{matrix}\frac{\bar{C_{L}} c_{0} r \rho^{2}}{4 \mu} s^{2} \sqrt{\mu^{2} + r^{2}} \left\lvert{\mu}\right\rvert \left(\dot{\eta}\right)^{2}\\0\\- \frac{\bar{C_{L}} c_{0}}{4} \rho^{2} s^{2} \sqrt{\mu^{2} + r^{2}} \left\lvert{\mu}\right\rvert \left(\dot{\eta}\right)^{2}\end{matrix}\right]$$

In [12]:

    

Dp_assump = Dp_expr.subs(assump_subs).to_matrix(frame_l).subs(assump_subs).subs(sub_adv).simplify()
Dp_assump


    

    
Out[12]:





$$\left[\begin{matrix}- \frac{\bar{C_{D}} c_{0}}{4} \rho^{2} s^{2} \sqrt{\mu^{2} + r^{2}} \left\lvert{\mu}\right\rvert \left(\dot{\eta}\right)^{2}\\0\\- \frac{\bar{C_{D}} c_{0} r \rho^{2}}{4 \mu} s^{2} \sqrt{\mu^{2} + r^{2}} \left\lvert{\mu}\right\rvert \left(\dot{\eta}\right)^{2}\end{matrix}\right]$$

In [13]:

    

L_l_x = sympy.integrate(Lp_assump[0], (r, 0, 1)).simplify()
L_l_x


    

    
Out[13]:





$$\frac{\bar{C_{L}} c_{0} \rho^{2} s^{2}}{12 \mu} \left(- \mu^{4} + \left(\mu^{2} + 1\right)^{\frac{3}{2}} \left\lvert{\mu}\right\rvert\right) \left(\dot{\eta}\right)^{2}$$

In [14]:

    

L_l_z = sympy.integrate(Lp_assump[2], (r, 0, 1))
L_l_z


    

    
Out[14]:





$$- \frac{\bar{C_{L}} c_{0}}{4} \rho^{2} s^{2} \left(\frac{\mu}{2} \sqrt{1 + \frac{1}{\mu^{2}}} + \frac{\mu \left\lvert{\mu}\right\rvert}{2} \operatorname{asinh}\left(\frac{1}{\left\lvert{\mu}\right\rvert}\right)\right) \left\lvert{\mu}\right\rvert \left(\dot{\eta}\right)^{2}$$

In [15]:

    

D_l_x = sympy.integrate(Dp_assump[0], (r, 0, 1)).simplify()
D_l_x


    

    
Out[15]:





$$- \frac{\bar{C_{D}} \mu}{8} c_{0} \rho^{2} s^{2} \left(\mu^{2} \operatorname{asinh}\left(\frac{1}{\left\lvert{\mu}\right\rvert}\right) + \sqrt{\mu^{2} + 1}\right) \left(\dot{\eta}\right)^{2}$$

In [16]:

    

D_l_z = sympy.integrate(Dp_assump[2], (r, 0, 1)).simplify()
D_l_z


    

    
Out[16]:





$$\frac{\bar{C_{D}} c_{0} \rho^{2} s^{2}}{12 \mu} \left(\mu^{4} - \left(\mu^{2} + 1\right)^{\frac{3}{2}} \left\lvert{\mu}\right\rvert\right) \left(\dot{\eta}\right)^{2}$$

In [17]:

    

F_A_l = sympy.Matrix([(L_l_x + D_l_x).simplify(), 0, (L_l_z + D_l_z).simplify()])
F_A_l


    

    
Out[17]:





$$\left[\begin{matrix}- \frac{c_{0} \rho^{2} s^{2} \left(\dot{\eta}\right)^{2}}{24 \mu} \left(3 \bar{C_{D}} \mu^{2} \left(\mu^{2} \operatorname{asinh}\left(\frac{1}{\left\lvert{\mu}\right\rvert}\right) + \sqrt{\mu^{2} + 1}\right) + 2 \bar{C_{L}} \left(\mu^{4} - \left(\mu^{2} + 1\right)^{\frac{3}{2}} \left\lvert{\mu}\right\rvert\right)\right)\\0\\\frac{c_{0} \rho^{2} s^{2} \left(\dot{\eta}\right)^{2}}{24 \mu} \left(2 \bar{C_{D}} \left(\mu^{4} - \left(\mu^{2} + 1\right)^{\frac{3}{2}} \left\lvert{\mu}\right\rvert\right) - 3 \bar{C_{L}} \mu^{2} \left(\mu^{2} \operatorname{asinh}\left(\frac{1}{\left\lvert{\mu}\right\rvert}\right) + \sqrt{\mu^{2} + 1}\right)\right)\end{matrix}\right]$$

In [18]:

    

const = {C_L_bar: 1.2, C_D_bar:0.5, rho:1.225, s:0.6, eta(t).diff(t):5, c_0:0.3}


    

In [19]:

    

N_to_g = 1000/9.8
two_wings = 2
sympy.plot(two_wings*N_to_g*(F_A_l[0]).subs(const),
           (mu,0,2), title='Thrust vs. $\mu$', xlabel='$\mu$', ylabel='Thrust, gm');


    

In [20]:

    

sympy.plot(-two_wings*N_to_g*(F_A_l[2]).subs(const),
           (mu,0,2), title='Lift vs. $\mu$', xlabel='$\mu$', ylabel='Lift, gm');


    

In [21]:

    

F_A_l.subs(const)[0]


    

    
Out[21]:





$$- \frac{1}{\mu} \left(0.40516875 \mu^{4} + 0.25323046875 \mu^{2} \left(\mu^{2} \operatorname{asinh}\left(\frac{1}{\left\lvert{\mu}\right\rvert}\right) + \sqrt{\mu^{2} + 1}\right) - 0.40516875 \left(\mu^{2} + 1\right)^{\frac{3}{2}} \left\lvert{\mu}\right\rvert\right)$$

In [22]:

    

f_thrust_mu = sympy.lambdify(mu, F_A_l.subs(const)[0], default_array=True)
f_lift_mu = sympy.lambdify(mu, -two_wings*N_to_g*F_A_l.subs(const)[2], default_array=True)
mu_max = scipy.optimize.newton(f_thrust_mu, 0.7)
lift_max = f_lift_mu(mu_max)
print('max lift: {:0.2f} gm @ mu = {:0.2f} -> V_T = {:0.2f} m/s'.format(
    lift_max, mu_max, mu_max*const[s]*const[eta(t).diff(t)]))


    

    




max lift: 501.13 gm @ mu = 1.25 -> V_T = 3.76 m/s